/*
 *  vector.cpp
 *  T3nsors2
 *
 *  Created by Michael Barriault on 10-06-16.
 *  Copyright 2010 University of Guelph. All rights reserved.
 *
 */

#include "vector.h"
#include "domain.h"
#include "scalar.h"
#include "others.h"
#include "macros.h"
#include <iostream>
#include <cstdlib>
#include <cmath>

namespace t3 {
	
	vector::vector(void) {
		return;
	}
	
	vector::~vector(void) {
		FOR(n,N) T[n]->release();
		delete[] T;
		O->release();
		if ( id != "" )
			print((string)"Deleting " + id);
	}
	
	vector* vector::init(domain* O, int N, string id) {
		vector* V;
		try {
			V = new vector;
//			V->T = (t3::scalar**)malloc(sizeof(scalar*)*N);
			V->T = new scalar*[N];
			FOR(n,N) V->T[n] = scalar::init(O, id+itos(n));
		}
		catch (...) {
			Alloc_Error(id);
			throw;
		}
		V->O = O;
		V->O->retain();
		V->N = N;
		V->id = id;
		V->count = 1;
		print((string)"Making new" + id);
		return V;
	}
	
	vector* vector::copy(vector* v) {
		if ( this->O == v->O ) FOR(n,this->N) FOR(o,this->O->size()) this->c(n)->at(o) = v->c(n)->at(o);
		return this;
	}
	
	vector* vector::retain(void) {
		print(id + "Going from " + itos(count) + " to ",true);
		count++;
		print(itos(count));
		return this;
	}

	void vector::release(void) {
		print(id + "Going from " + itos(count) + " to ",true);
		count--;
		print(itos(count));
		if ( count <= 0 ) delete this;
	}
	
	scalar*& vector::c(int a) {
		return T[a];
	}
	
	vector* diff(scalar* x) {
		vector* dx = vector::init(x->O, x->O->N, (string)"d"+x->id);
		FOR(p,dx->N) {
			int dg = x->O->upto(p);
			int dc = x->O->upto(p-1);
			int G = x->O->size()/dg;
			int C = x->O->at(p);
			FOR(g,G) {
				FOR(m,dc) {
					int i = g*dg + m;
					dx->c(p)->at(i) = (x->at(i+dc)-x->at(i))/dx->O->d(p);
					//dx->c(p)->at(i) = (-3*x->at(i)+4*x->at(i+dc)-x->at(i+2*dc))/(2*dx->O->d(p));
					for ( int c=1; c<C-1; c++ ) {
						i = g*dg + m + c*dc;
						dx->c(p)->at(i) = (x->at(i+dc)-x->at(i-dc))/(2*dx->O->d(p));
					}
					i = g*dg + m + (C-1)*dc;
					dx->c(p)->at(i) = (x->at(i)-x->at(i-dc))/dx->O->d(p);
					//dx->c(p)->at(i) = (x->at(i-2*dc)-4*x->at(i-dc)+3*x->at(i))/(2*dx->O->d(p));
				}
			}
		}
		return dx;
	}

	scalar* lapfour(scalar* x) {
		vector* dx = vector::init(x->O, x->O->N, (string)"d4"+x->id);
		FOR(p,dx->N) {
			int dg = x->O->upto(p);
			int dc = x->O->upto(p-1);
			int G = x->O->size()/dg;
			int C = x->O->at(p);
			FOR(g,G) {
				FOR(m,dc) {
					int i = g*dg + m;
					dx->c(p)->at(i) = 0;
					dx->c(p)->at(i+dc) = 0;
					for ( int c=2; c<C-2; c++ ) {
						i = g*dg + m + c*dc;
						dx->c(p)->at(i) = (x->at(i-2*dc)-4*x->at(i-dc)+6*x->at(i)-4*x->at(i+dc)+x->at(i+2*dc))/pow(x->O->d(p),4.);
					}
					i = g*dg + m + (C-1)*dc;
					dx->c(p)->at(i-dc) = 0;
					dx->c(p)->at(i) = 0;
				}
			}
		}
		scalar *Dx = scalar::init(x->O, dx->id);
		FOR(a,dx->N) FOR(o,x->O->size()) Dx->at(o) += dx->c(a)->at(o);
		dx->release();
		return Dx;
	}
	
	vector* lapfour(vector* x) {
		vector* dx = vector::init(x->O, x->N, (string)"d4"+x->id);
		FOR(p,dx->N) {
			dx->c(p)->release();
			dx->c(p) = lapfour(x->c(p));
		}
		return dx;
	}
	
	scalar* dotProduct(vector* x, vector* y) {
		scalar* xy = scalar::init(x->O, x->id+"*"+y->id);
		if ( x->O == y->O and x->N == y->N )
			FOR(a,x->N) FOR(o,x->O->size()) xy->at(o) += x->c(a)->at(o)*y->c(a)->at(o);
		return xy;
	}

	vector* linC(double a, vector* x, double b, vector* y, string id) {
		vector* z = vector::init(x->O, x->N, id);
		if ( x->O == y->O and x->N == y->N )
			FOR(a,x->N) FOR(o,x->O->size()) z->c(a)->at(o) = a*x->c(a)->at(o) + b*y->c(a)->at(o);
		return z;
	}

	vector* linC(double a, vector* x, double b, vector* y, double c, vector* z, string id) {
		vector* w = vector::init(x->O, x->N, id);
		if ( x->O == y->O and x->N == y->N and y->O == z->O and y->N == z->N )
			FOR(a,x->N) FOR(o,x->O->size()) w->c(a)->at(o) = a*x->c(a)->at(o) + b*y->c(a)->at(o) + c*z->c(a)->at(o);
		return w;
	}
}
